Functional analysis - bacteria (ggpicrust2)

Here is an attempt at functional analysis using R

QZA files generated in Picrust2 were imported using qiime2r and converted to data frames

Non-significant KO abundance analysis (ggpicrust2)

Several ggpicrust2 tests failed due to no significant pathways or biomarkers. The code for those is in the tabs below. These will be removed for the final document

Treatment - aggregated

#ko_abund.res <- ggpicrust2(data = ko_abund.df,
#                           metadata = metadata.df,
#                           group="Treatment",
#                           pathway="KO",
#                           ko_to_kegg = TRUE,
#                           p.adjust="fdr",
#                           order="pathway_class",
#                           p_values_bar=TRUE,
#                           x_lab="pathway_name")

# Failed - no significant biomarkers identified

Treatment - unaggregated

#ko_abund.res.nonagg <- ggpicrust2(data = ko_abund.df.nonagg,
#                           metadata = metadata.df.nonagg,
#                           group="Treatment",
#                           pathway="KO",
#                           ko_to_kegg = TRUE,
#                           p.adjust="fdr",
#                           order="pathway_class",
#                           p_values_bar=TRUE,
#                           x_lab="pathway_name")

# Failed - no significant biomarkers identified

Treatment - site by site

# Harley Farms
#ko_abund.res.HF <- ggpicrust2(data = ko_abund.HF,
#                           metadata = metadata_site.split$HARLEY_FARMS_new2b,
#                           group="Treatment",
#                           pathway="KO",
#                           ko_to_kegg = TRUE,
#                           p.adjust="fdr",
#                           order="pathway_class",
#                           p_values_bar=TRUE,
#                           x_lab="pathway_name")

# Windmill Farm
#ko_abund.res.WF <- ggpicrust2(data = ko_abund.WF,
#                           metadata = metadata_site.split$WINDMILL_farm,
#                           group="Treatment",
#                           pathway="KO",
#                           ko_to_kegg = TRUE,
#                           p.adjust="fdr",
#                           order="pathway_class",
#                           p_values_bar=TRUE,
#                           x_lab="pathway_name")

# Castle Field West
#ko_abund.res.CW <- ggpicrust2(data = ko_abund.CW,
#                           metadata = metadata_site.split$Castle_Field_West_W,
#                           group="Treatment",
#                           pathway="KO",
#                           ko_to_kegg = TRUE,
#                           p.adjust="fdr",
#                           order="pathway_class",
#                           p_values_bar=TRUE,
#                           x_lab="pathway_name")

# Jemma 6
#ko_abund.res.J6 <- ggpicrust2(data = ko_abund.J6,
#                           metadata = metadata_site.split$Jemma_6,
#                           group="Treatment",
#                           pathway="KO",
#                           ko_to_kegg = TRUE,
#                           p.adjust="fdr",
#                           order="pathway_class",
#                           p_values_bar=TRUE,
#                           x_lab="pathway_name")

# Jemma 9
#ko_abund.res.J9 <- ggpicrust2(data = ko_abund.J9,
#                           metadata = metadata_site.split$Jemma_9,
#                           group="Treatment",
#                           pathway="KO",
#                           ko_to_kegg = TRUE,
#                           p.adjust="fdr",
#                           order="pathway_class",
#                           p_values_bar=TRUE,
#                           x_lab="pathway_name")

# Lardon Chase
#ko_abund.res.LC <- ggpicrust2(data = ko_abund.LC,
#                           metadata = metadata_site.split$Lardon_Chase,
#                           group="Treatment",
#                           pathway="KO",
#                           ko_to_kegg = TRUE,
#                           p.adjust="fdr",
#                           order="pathway_class",
#                           p_values_bar=TRUE,
#                           x_lab="pathway_name")

# No statistically significant pathways in any of the sites when comparing control and shelter samples

Site comparisons - shelter samples

#ko_abund.res.shelt <- ggpicrust2(data = ko_abund.shelt,
#                           metadata = metadata_treat.split$Drought,
#                           group="Site",
#                           pathway="KO",
#                           ko_to_kegg = TRUE,
#                           p.adjust="fdr",
#                           order="pathway_class",
#                           p_values_bar=TRUE,
#                           x_lab="pathway_name")

# Failed - no statistically significant biomarkers in the dataset

Significant pathways from KO abundance analysis (ggpicrust2) - control samples

Error barplot - Kruskal-Wallis

Error barplot - GLM (top 20)

Error barplot - GLM (top 50)

PCA plot - KO abundance

Heatmap - Kegg pathway expression in samples relative to mean abundances

GLM

## The Sample Names in order from left to right are:
## P2C1, P2C2, P2C3, P19C1, P19C2, P19C3, P54C1, P54C2, P54C3, P55C1, P55C2, P55C3, P67C1, P67C2, P67C3, P20C1, P20C2, P20C3
## The Group Levels in order from left to right are:
## CW, CW, CW, HF, HF, HF, J6, J6, J6, J9, J9, J9, LC, LC, LC, WF, WF, WF

Kruskal-Wallis

## The Sample Names in order from left to right are:
## P2C1, P2C2, P2C3, P19C1, P19C2, P19C3, P54C1, P54C2, P54C3, P55C1, P55C2, P55C3, P67C1, P67C2, P67C3, P20C1, P20C2, P20C3
## The Group Levels in order from left to right are:
## CW, CW, CW, HF, HF, HF, J6, J6, J6, J9, J9, J9, LC, LC, LC, WF, WF, WF

Most abundant KEGG pathways in control samples

Non-significant MetaCyc pathway abundance analysis (ggpicrust2)

Several ggpicrust2 tests failed due to no significant pathways. The code for those is in the tabs below. These will be removed for the final document

Treatment - aggregated

#path_abund.res <- ggpicrust2(data = path_abund.df,
#                             metadata = metadata.df,
#                             group="Treatment",
#                             pathway="MetaCyc",
#                             ko_to_kegg = FALSE,
#                             p.adjust="fdr",
#                             order="group",
#                             p_values_bar = TRUE,
#                             x_lab="description")

# Failed - no differentially abundant pathways

Treatment - unaggregated

#path_abund.res.nonagg <- ggpicrust2(data = path_abund.df.nonagg,
#                             metadata = metadata.df.nonagg,
#                             group="Treatment",
#                             pathway="MetaCyc",
#                             ko_to_kegg = FALSE,
#                             p.adjust="fdr",
#                             order="group",
#                             p_values_bar = TRUE,
#                             x_lab="description")

# Failed - no significantly differentially abundant pathways

Treatment - site by site

# Harley Farms
#path_abund.res.HF <- ggpicrust2(data = path_abund.HF,
#                           metadata = metadata_site.split$HARLEY_FARMS_new2b,
#                           group="Treatment",
#                           pathway="MetaCyc",
#                           ko_to_kegg = FALSE,
#                           p.adjust="fdr",
#                           order="group",
#                           p_values_bar=TRUE,
#                           x_lab="description")

# Windmill Farm
#path_abund.res.WF <- ggpicrust2(data = path_abund.WF,
#                           metadata = metadata_site.split$WINDMILL_farm,
#                           group="Treatment",
#                           pathway="MetaCyc",
#                           ko_to_kegg = FALSE,
#                           p.adjust="fdr",
#                           order="group",
#                           p_values_bar=TRUE,
#                           x_lab="description")

# Castle Field West
#path_abund.res.CW <- ggpicrust2(data = path_abund.CW,
#                           metadata = metadata_site.split$Castle_Field_West_W,
#                           group="Treatment",
#                           pathway="MetaCyc",
#                           ko_to_kegg = FALSE,
#                           p.adjust="fdr",
#                           order="group",
#                           p_values_bar=TRUE,
#                           x_lab="description")

# Jemma 6
#path_abund.res.J6 <- ggpicrust2(data = path_abund.J6,
#                           metadata = metadata_site.split$Jemma_6,
#                           group="Treatment",
#                           pathway="MetaCyc",
#                           ko_to_kegg = FALSE,
#                           p.adjust="fdr",
#                           order="group",
#                           p_values_bar=TRUE,
#                           x_lab="description")

# Jemma 9
#path_abund.res.J9 <- ggpicrust2(data = path_abund.J9,
#                           metadata = metadata_site.split$Jemma_9,
#                           group="Treatment",
#                           pathway="MetaCyc",
#                           ko_to_kegg = FALSE,
#                           p.adjust="fdr",
#                           order="group",
#                           p_values_bar=TRUE,
#                           x_lab="description")

# Lardon Chase
#path_abund.res.LC <- ggpicrust2(data = path_abund.LC,
#                           metadata = metadata_site.split$Lardon_Chase,
#                           group="Treatment",
#                           pathway="MetaCyc",
#                           ko_to_kegg = FALSE,
#                           p.adjust="fdr",
#                           order="group",
#                           p_values_bar=TRUE,
#                           x_lab="description")

# No statistically significant MetaCyc pathways in any of the sites when comparing control and shelter samples

Site comparisons - control samples

#path_abund.res.ctrl <- ggpicrust2(data = path_abund.ctrl,
#                           metadata = metadata_treat.split$Control,
#                           group="Site",
#                           pathway="MetaCyc",
#                           ko_to_kegg = FALSE,
#                           p.adjust="fdr",
#                           order="group",
#                           p_values_bar=TRUE,
#                           x_lab="description")

# Failed due to lack of statistically significant features

Most abundant MetaCyc pathways in shelter samples

PCA plot - MetaCyc pathway

Significant MetaCyc pathways from pathway abundance analysis (ggpicrust2) - shelter samples

Error barplot - Kruskal-Wallis

Error barplot - GLM (top 3)

Heatmap - Metacyc pathway expression in samples relative to mean abundances

GLM

## The Sample Names in order from left to right are:
## P2S1, P2S2, P2S3, P19S1, P19S2, P19S3, P54S1, P54S2, P54S3, P55S1, P55S2, P55S3, P67S1, P67S2, P67S3, P20S1, P20S2, P20S3
## The Group Levels in order from left to right are:
## CW, CW, CW, HF, HF, HF, J6, J6, J6, J9, J9, J9, LC, LC, LC, WF, WF, WF

Kruskal-Wallis

## The Sample Names in order from left to right are:
## P2S1, P2S2, P2S3, P19S1, P19S2, P19S3, P54S1, P54S2, P54S3, P55S1, P55S2, P55S3, P67S1, P67S2, P67S3, P20S1, P20S2, P20S3
## The Group Levels in order from left to right are:
## CW, CW, CW, HF, HF, HF, J6, J6, J6, J9, J9, J9, LC, LC, LC, WF, WF, WF

Functional analysis - fungi (FUNGuildR)

Top 10 most abundant guilds

Control

Shelter

Guilds grouped by site-treatment

Guilds grouped by treatment